###
# Replication for Appendix: Simmons, Lloyd and Stewart
###
# Code runs subject to the assumption your working directory is
# set to the location of this script.

# A Note on the Data: 
# The data was collected to run through criminalizations in 2008
# we use a particular hack of listing the criminalization date
# of countries that we know have not yet criminalized as 2012 and
# list as NA those whose criminalization status is unknown entirely.
# This works for this code but should obviously be addressed if
# one is attempting to update the data.

# This appendix starts with a dataset which is used to do the
# modeling in the main replication code.  You can grab it there.

library(survival)
library(foreign)
library(data.table)
library(Greg)

# A function to make tables we can use in appendix or for the paper
# simply calculates hazard ratios of a cox model along with CI's and 
# p-values.  Either returns the table or the xtable 
returnTable <- function(mod,rownames=NULL,xtable=FALSE, digits=3, caption=NULL) {
  tab <- summary(mod)$coefficients[,c(1,2,6)]
  Lower <- exp(tab[,1] + qnorm(.025)*tab[,2])
  Upper <- exp(tab[,1] + qnorm(.975)*tab[,2])
  HR <- exp(tab[,1])
  tab <- data.frame(HR, Lower, Upper, pvalue=tab[,3])
  rownames(tab) <- rownames
  colnames(tab) <- c("Hazard Ratio", "Lower 95% CI", "Upper 95% CI", "p-value")
  if(xtable) {
    caption <- sprintf("%s \\newline %s", caption, diagnose(mod, xtable=TRUE))
    return(xtable::xtable(tab, digits=digits, caption=caption))
  } else {
    return(round(tab, digits))
  }
}

#prints out number of observations, number of events,
#and the cox non-proportional hazards tests of Grambsch and Therneau (1994)
#for a cox model
diagnose <- function(mod, xtable=FALSE) {
  n <- summary(mod)$n
  ne <- summary(mod)$nevent
  p.nph <- cox.zph(mod)$table[,3] #grab p-values
  global <- p.nph[length(p.nph)] < .05
  nph <- any(p.nph < .05)
  if(xtable) {
   #print text to be used with xtable
    sprintf("%i complete cases (%i events). There were %i p-values less than .05 in the NPH tests with a global p-valiue of %.3f",
            n, ne, sum(p.nph < .05), p.nph[length(p.nph)])
   #print(xtable::xtable(matrix(c(n,ne,sum(p.nph < .05)), ncol=1, 
   #                       dimnames=list(c("N Complete Cases", "Events", "N NPH Prop p<.05"))),
  #                digits=0),include.colnames=FALSE)
  } else {
    cat(sprintf("N Complete Cases: \t %i \nN Events: \t \t %i \nAny NPH Prop p<.05: \t %s \nGlobal NPH Prop p<.05: \t %s",
              n, ne, nph, global))
  }
}

#Make a little function that will add p-values
#to an ordinal probit table created by MASS.  Return
#a rounded table out to 3 digits.  Also provides number of obs
op_summary <- function(mod) {
  m.coef <- data.frame(coef(summary(mod)))
  m.coef$pval <- round((pnorm(abs(m.coef$t.value), lower.tail = FALSE) * 2),2)
  list(coef=round(m.coef,3), n=mod$n)
}

make_na_data <- function(form, data) {
  mini <- data[,all.vars(as.formula(form)),with=FALSE]
  na.omit(mini)
}

############
# Generate Tables/Figures
############
load("HTreadytomodel.RData")

#########
# Comparison asinh comparison figure

form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2) + uspressure_1 + pspline(rule_of_law, df=3) + rattocprot +
  cluster(country)
t1m1 <- coxph(form, data=data, x=TRUE)
form <- Surv(t0,t,fail) ~ asinh(Crimsum) + uspressure_1 + pspline(rule_of_law, df=3) + rattocprot +
  cluster(country)
t1m1a <- coxph(form, data=data, x=TRUE)

pdf(file = "figures/CrimFunctionCompare.pdf", width = 8, height = 8, family = "Helvetica") 
par(mfrow=c(2,2))
Greg::plotHR(t1m1, xlim=c(0,10), ylim=c(.5, 4), rug="ticks", main="Spatial Spline (Zoomed)", xlab="Number of Roads to Criminalized Neighbors") 
Greg::plotHR(t1m1, xlim=c(0,60), ylim=c(.5, 4), rug="ticks", main="Spatial Spline", xlab="Number of Roads to Criminalized Neighbors")
Greg::plotHR(t1m1a, xlim=c(0,10), ylim=c(.5, 4), rug="ticks", main="Spatial Inverse Hyperbolic Sine", xlab="Number of Roads to Criminalized Neighbors") 
Greg::plotHR(t1m1a, xlim=c(0,60), ylim=c(.5, 4), rug="ticks", main="Spatial Inverse Hyperbolic Sine", xlab="Number of Roads to Criminalized Neighbors")
dev.off()



#########
##Table 1, Model 1
form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2) + uspressure_1 + 
  pspline(rule_of_law, df=3) + rattocprot + cluster(country)
t1m1 <- coxph(form, data=data, x=TRUE)

sink(file = "tex/T1M1.tex")
returnTable(t1m1, xtable=TRUE,
                   rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                              "US Pressure", "Rule of Law",
                              "Rule of Law Nonlinear", "Ratification"),
                   caption="Table 1, Model 1.")
sink()

pdf(file = "figures/t1m1.pdf", width =10, height = 4, family = "Helvetica") 
par(mfrow=c(1,2))
Greg::plotHR(t1m1, xlim=c(0,60), ylim=c(.5, 4), rug="ticks", main="Spatial Spline", xlab="Number of Roads to Criminalized Neighbors")
Greg::plotHR(t1m1,term=3, rug="ticks", xlim=c(-2.5,2.5), xlab="Rule of Law",
             main="Rule of Law")
dev.off()

#########
#Table 1, Model 2
form <- Surv(t0,t,fail) ~ pspline(Crimsum, df=2) + uspressure_1 + pspline(rule_of_law, df=3) + 
  rattocprot + pspline(USaid_gdp, df=2) +
  useimfcr + 
  pspline(us_trade_share_cap, df=2) + 
  pspline(eu_trade_share_cap, df=2) +
  cluster(country)
t1m2 <- coxph(form, data=data, x=TRUE)

sink(file = "tex/T1M2.tex")
returnTable(t1m2, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                             "US Pressure", "Rule of Law",
                             "Rule of Law Nonlinear", "Ratification",
                             "US aid GDP", "US Aid Nonlinear",
                             "IMF Credit", "US Trade Share", 
                             "US Trade Share Nonlinear",
                             "EU Trade Share", "EU Trade Share Nonlinear"),
            xtable=TRUE, caption="Table 1, Model 2.")
sink()

pdf(file="figures/T1M2.pdf", width=15, height=8)
par(mfrow=c(2,3)) 
termplot(t1m2, term=1, main="Spatial Spline", xlab="Number of Roads to Criminalized Neighbors", se=TRUE)
termplot(t1m2, term=3, main="Rule of Law", xlab="Rule of Law", se=TRUE)
termplot(t1m2, term=5, main="US Aid", xlab="US Aid by GDP", se=TRUE)
termplot(t1m2, term=7, main="US Trade Share", xlab="US Trade Share", se=TRUE)
termplot(t1m2, term=8, main="EU Trade Share", xlab="EU Trade Share", se=TRUE)
dev.off()


#######
#Table 1, Model 3

form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2) + uspressure_1 + 
  pspline(rule_of_law, df=3) + rattocprot + 
  pspline(ichild_labor, df=3) + mid + 
  pspline(remittances,df=3) +
  cluster(country)


t1m3 <- coxph(form, data=data, x=TRUE)

sink(file = "tex/T1M3.tex")
tab <- returnTable(t1m3, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                                    "US Pressure", "Rule of Law",
                                    "Rule of Law Nonlinear", "Ratification",
                                    "Child Labor", "Child Labor Nonlinear",
                                    "Middle Income Country",
                                    "Remittances", "Remittances Nonlinear"),
                   xtable=TRUE, caption="Table 1, Model 3.")
tab
sink()

pdf(file = "figures/T1M3.pdf", width = 10, height = 8, family = "Helvetica") 
par(mfrow=c(2,2))
termplot(t1m3, term=1, se=TRUE, main="Spatial Spline", xlab="Number of Roads to Criminalized Neighbors")
termplot(t1m3, term=3, se=TRUE, xlab="Rule of Law", main="Rule of Law")
termplot(t1m3, term=5, se=TRUE, xlab="Child Labor", main="Child Labor")
termplot(t1m3, term=7, se=TRUE, xlab="Remittances", main="Remittances")
dev.off()

#########
#Table 1, Model 4
form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2) + uspressure_1 + 
  pspline(rule_of_law, df=3) + rattocprot + 
  pspline(ichild_labor,df=3) + mid +
  pspline(par_women_lower_, df=3) +
  cluster(country)
t1m4 <- coxph(form, data=data, x=TRUE)

sink(file = "tex/T1M4.tex")
tab <- returnTable(t1m4, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                                    "US Pressure", "Rule of Law",
                                    "Rule of Law Nonlinear", "Ratification",
                                    "Child Labor", "Child Labor Nonlinear",
                                    "Middle Income Country",
                                    "Share Women in Parliament",
                                    "Share Women Nonlinear"),
                   xtable=TRUE, caption="Table 1, Model 4.")
tab
sink()

pdf(file = "figures/T1M4.pdf", width = 10, height = 8, family = "Helvetica") 
par(mfrow=c(2,2))
termplot(t1m4, term=1, se=TRUE, main="Spatial Spline", xlab="Number of Roads to Criminalized Neighbors")
termplot(t1m4, term=3, se=TRUE, xlab="Rule of Law", main="Rule of Law")
termplot(t1m4, term=5, se=TRUE, xlab="Child Labor", main="Child Labor")
termplot(t1m4, term=7, se=TRUE, xlab="Share of Women in Parliament", main="Women's Participation")
dev.off()



#########
#Table 1, Model 5
form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2) + uspressure_1 + 
  pspline(rule_of_law, df=3) + rattocprot+ 
  pspline(latentmean, df=3) +
  cluster(country)

t1m5 <- coxph(form, data=data, x=TRUE)

sink(file = "tex/T1M5.tex")
tab <- returnTable(t1m5, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                                    "US Pressure", "Rule of Law",
                                    "Rule of Law Nonlinear", "Ratification",
                                    "Respect For Human Rights",
                                    "Respect for HR Nonlinear"),
                   xtable=TRUE, caption="Table 1, Model 5.")
tab
sink()

pdf(file = "figures/T1M5.pdf", width = 10, height = 8, family = "Helvetica") 
par(mfrow=c(2,2))
termplot(t1m5, term=1, se=TRUE, main="Spatial Spline", xlab="Number of Roads to Criminalized Neighbors")
termplot(t1m5, term=3, se=TRUE, xlab="Rule of Law", main="Rule of Law")
termplot(t1m5, term=5, se=TRUE, xlab="Respect for Human Rights", main="Respect for HR Score")
dev.off()

############
#Table 1, Model 6
form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2) + uspressure_1 + 
  pspline(rule_of_law, df=3) + rattocprot + 
  pspline(crimsumany_HRIO, df=3) +
  cluster(country)
t1m6 <- coxph(form, data=data, x=TRUE)

sink(file = "tex/T1M6.tex")
tab <- returnTable(t1m6, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                                    "US Pressure", "Rule of Law",
                                    "Rule of Law Nonlinear", "Ratification",
                                    "Shared HROrg Crim",
                                    "Shared HROrg Crim Nonlinear"),
                   xtable=TRUE, caption= "Table 1, Model 6.")
tab
sink()

pdf(file = "figures/T1M6.pdf", width = 10, height = 8) 
par(mfrow=c(2,2))
termplot(t1m6, term=1, se=TRUE, main="Spatial Spline", xlab="Number of Roads to Criminalized Neighbors")
termplot(t1m6, term=3, se=TRUE, xlab="Rule of Law", main="Rule of Law")
termplot(t1m6, term=5, se=TRUE, xlab="Shared HROrg", main="Shared HR Organizations")
dev.off()

#code for the p-values table referenced in the appendix
cox.zph(t1m6)

##########
# Table 2: Robustness of Diffusion via Roads: Other "Proximity" Measures
##########

#Model 1: Policy Weighted by Trade Partners
form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2)+ uspressure_1 + 
  pspline(rule_of_law, df=3) + pspline(crim_ht_w1, df=2) +
  cluster(country)
t2m1 <- coxph(form, data=data, x=TRUE)

sink(file = "tex/T2M1.tex")
tab <- returnTable(t2m1, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                                    "US Pressure", "Rule of Law",
                                    "Rule of Law Nonlinear", 
                                    "Policy Weighted by Trade Partners",
                                    "Policy Weighted by Trade Nonlinear"),
                   xtable=TRUE, caption= "Table 2, Model 1.")
tab
sink()

pdf(file = "figures/T2M1.pdf", width = 10, height = 8) 
par(mfrow=c(2,2))
termplot(t2m1, term=1, se=TRUE, main="Spatial Spline", xlab="Number of Roads to Criminalized Neighbors")
termplot(t2m1, term=3, se=TRUE, xlab="Rule of Law", main="Rule of Law")
termplot(t2m1, term=4, se=TRUE, xlab="Policy Weighted by Trade Partners", main="Trade")
dev.off()

#Model 2: Regional Press Stories on HT
form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2)+ uspressure_1 + 
  pspline(rule_of_law, df=3) + pspline(laglogreg_arts2, df=2) +
  cluster(country)
t2m2 <- coxph(form, data=data, x=TRUE)


sink(file = "tex/T2M2.tex")
tab <- returnTable(t2m2, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                                    "US Pressure", "Rule of Law",
                                    "Rule of Law Nonlinear", 
                                    "Regional News",
                                    "Regional News Nonlinear"),
                   xtable=TRUE, caption= "Table 2, Model 2.")
tab
sink()

pdf(file = "figures/T2M2.pdf", width = 10, height = 8) 
par(mfrow=c(2,2))
termplot(t2m2, term=1, se=TRUE, main="Spatial Spline", xlab="Number of Roads to Criminalized Neighbors")
termplot(t2m2, term=3, se=TRUE, xlab="Rule of Law", main="Rule of Law")
termplot(t2m2, term=4, se=TRUE, xlab="Regional News", main="Regional News")
dev.off()

#Model 3: Developmental Emulation
form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2)+ uspressure_1 + 
  pspline(rule_of_law, df=3) + pspline(devlevDens, df=2) + 
  cluster(country)
t2m3 <- coxph(form, data=data, x=TRUE)

sink(file = "tex/T2M3.tex")
tab <- returnTable(t2m3, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                                    "US Pressure", "Rule of Law",
                                    "Rule of Law Nonlinear", 
                                    "Density in DevLevel",
                                    "Density in DevLevel Nonlinear"),
                   xtable=TRUE, caption= "Table 2, Model 3.")
tab
sink()

pdf(file = "figures/T2M3.pdf", width = 10, height = 8) 
par(mfrow=c(2,2))
termplot(t2m3, term=1, se=TRUE, main="Spatial Spline", xlab="Number of Roads to Criminalized Neighbors")
termplot(t2m3, term=3, se=TRUE, xlab="Rule of Law", main="Rule of Law")
termplot(t2m3, term=4, se=TRUE, xlab="Density in Development Level", main="Development")
dev.off()

#Model 4: Civilizational Emulation
form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2)+ uspressure_1 + 
  pspline(rule_of_law, df=3) + pspline(civDens, df=2)  + cluster(country)
t2m4 <- coxph(form, data=data, x=TRUE)


sink(file = "tex/T2M4.tex")
tab <- returnTable(t2m4, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                                    "US Pressure", "Rule of Law",
                                    "Rule of Law Nonlinear",
                                    "Civilizational Crim",
                                    "CivCrim Nonlinear"),
                   xtable=TRUE, caption= "Table 2, Model 4.")
tab
sink()

pdf(file = "figures/T2M4.pdf", width = 10, height = 8) 
par(mfrow=c(2,2))
Greg::plotHR(t2m4, xlim=c(0,60), ylim=c(.5, 4), rug="ticks", main="Spatial Spline", xlab="Number of Roads to Criminalized Neighbors")
Greg::plotHR(t2m4,term=3, rug="ticks", xlim=c(-2.5,2.5), xlab="Rule of Law",
             main="Rule of Law")
Greg::plotHR(t2m4,term=4, rug="ticks", xlim=c(0,60), xlab="Civilizational Neighbor Criminalizations",
             main="Civilizational Emulation")
dev.off()

####
#Model 5: Legal Families
####

#NB: there was a mistake made in early processing where legal density
#    was calculated and then multiplied by 100 (to make it a percentage)
#    however this was done twice. Dividing back by a 100 gets us to nice
#    numbers again.  We adjust here.
data$legalDens <- data$legalDens/100

form <- Surv(t0,t,fail ) ~ pspline(Crimsum,df=2)+ uspressure_1 + 
  pspline(rule_of_law, df=3) + pspline(legalDens, df=2)  + cluster(country)
t2m5 <- coxph(form, data=data, x=TRUE)

sink(file = "tex/T2M5.tex")
tab <- returnTable(t2m5, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                                    "US Pressure", "Rule of Law",
                                    "Rule of Law Nonlinear",
                                    "Legal Crim",
                                    "LegCrim Nonlinear"),
                   xtable=TRUE, caption= "Table 2, Model 5.")
tab
sink()


pdf(file = "figures/T2M5.pdf", width = 10, height = 8, family = "Helvetica") 
par(mfrow=c(2,2))
Greg::plotHR(t2m5, xlim=c(0,60), ylim=c(.5, 4), rug="ticks", main="Spatial Spline", xlab="Number of Roads to Criminalized Neighbors")
Greg::plotHR(t2m5,term=3, rug="ticks", xlim=c(-2.5,2.5), xlab="Rule of Law",
             main="Rule of Law")
Greg::plotHR(t2m5,term=4, rug="ticks", xlim=c(0,50), xlab="Legal Neighbor Criminalizations",
             main="Legal Emulation")
dev.off()

#########
# Table 3 (Model 2) / Figure 5: HT and ML Comparison
#########
##
#Human Trafficking DV
##
form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2)+ 
  pspline(rule_of_law, df=3) + devlev  + cluster(country)
t3m2a <- coxph(form, data=data, x=TRUE)

sink(file = "tex/T3M2a.tex")
tab <- returnTable(t3m2a, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                                     "Rule of Law",
                                     "Rule of Law Nonlinear",
                                     "DevLev"),
                   xtable=TRUE, caption= "Table 3, Model 2a.")
tab
sink()

pdf(file = "figures/T3M2a.pdf", width = 10, height = 4, family = "Helvetica") 
par(mfrow=c(1,2))
Greg::plotHR(t3m2a, xlim=c(0,60), ylim=c(.5, 4), rug="ticks", main="Spatial Spline", xlab="Number of Roads to Criminalized Neighbors")
Greg::plotHR(t3m2a,term=2, rug="ticks", xlim=c(-2.5,2.5), xlab="Rule of Law",
             main="Rule of Law")
dev.off()

##
#Money laundering 
##
# we now need to load a new dataset and recreate
# the survival model variables

htdata <- data
load("CompleteDataWithWeights.RData")
data <- data.table(data)

#Start timer at 1990
data$t <- data$year - 1990
data$t0 <- data$year - 1991
subset <- data[year>=1990,]
setkey(subset, country, year)

#generate lags and fail conditions
subset[,crim_ml_lag1 := shift(ml_crim_lax), by=country]
subset[,fail := ifelse(crim_ml_lag1==0 & ml_crim_lax==1, 1, 0)]
subset[,drop := ifelse(crim_ml_lag1==1, 1, 0)]
#drop out those who have already failed
subset <- subset[drop!=1,]
subset$drop <- NULL
#rename
MLdata <- subset
data <- htdata

form <- Surv(t0,t,fail) ~ pspline(MLsum,df=2)+ 
  pspline(rule_of_law, df=3)+ devlev + cluster(country)
t3m2b <- coxph(form, data=MLdata, x=TRUE)


sink(file = "tex/T3M2b.tex")
tab <- returnTable(t3m2b, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                                     "Rule of Law",
                                     "Rule of Law Nonlinear",
                                     "DevLev"),
                   xtable=TRUE, caption= "Table 3, Model 2b.")
tab
sink()


pdf(file = "figures/T3M2b.pdf", width = 10, height = 4, family = "Helvetica") 
par(mfrow=c(1,2))
Greg::plotHR(t3m2b, xlim=c(0,60), ylim=c(.5, 12), rug="ticks", 
             main="Money Laundering", 
             xlab="Number of Roads to Criminalized Neighbors")
Greg::plotHR(t3m2b,term=2, rug="ticks", xlim=c(-2.5,2.5), xlab="Rule of Law",
             main="Rule of Law")
dev.off()

#code for the p-values table referenced in the appendix
cox.zph(t3m2b)


########
# Table 4/Figure 6: By Exposure Category
########

#now we are back to human trafficking so we load back up
# that data
load("HTreadytomodel.RData")

##
#Destination Countries
##
form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2)+ uspressure_1 + 
  pspline(rule_of_law, df=3) + rattocprot + cluster(country)
t4m1 <- coxph(form, data=data[data$destination==1,], x=TRUE)

sink(file = "tex/T4M1.tex")
tab <- returnTable(t4m1, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                                    "US Pressure", "Rule of Law",
                                    "Rule of Law Nonlinear",
                                    "Ratification"),
                   xtable=TRUE, caption= "Table 4, Model 1.")
tab
sink()

pdf(file = "figures/T4M1.pdf", width = 10, height = 4, family = "Helvetica") 
par(mfrow=c(1,2))
Greg::plotHR(t4m1, xlim=c(0,60), ylim=c(.5, 4), rug="ticks", main="Spatial Spline", xlab="Number of Roads to Criminalized Neighbors")
Greg::plotHR(t4m1,term=3, rug="ticks", xlim=c(-2.5,2.5), xlab="Rule of Law",
             main="Rule of Law")
dev.off()

##
#Origin Countries
##
form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2)+ uspressure_1 + 
  pspline(rule_of_law, df=3) + rattocprot + cluster(country)
t4m2 <- coxph(form, data=data[data$origin==1,], x=TRUE)


sink(file = "tex/T4M2.tex")
tab <- returnTable(t4m2, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                                    "US Pressure", "Rule of Law",
                                    "Rule of Law Nonlinear",
                                    "Ratification"),
                   xtable=TRUE, caption= "Table 4, Model 1.")
tab
sink()

pdf(file = "figures/T4M2.pdf", width = 10, height = 4, family = "Helvetica") 
par(mfrow=c(1,2))
Greg::plotHR(t4m2, xlim=c(0,60), ylim=c(.5, 4), rug="ticks", main="Spatial Spline", xlab="Number of Roads to Criminalized Neighbors")
Greg::plotHR(t4m2,term=3, rug="ticks", xlim=c(-2.5,2.5), xlab="Rule of Law",
             main="Rule of Law")
dev.off()


##
#Transit Countries
##
form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2)+ uspressure_1 + 
  pspline(rule_of_law, df=3) + rattocprot + cluster(country)
t4m3 <- coxph(form, data=data[data$transit==1,], x=TRUE)


sink(file = "tex/T4M3.tex")
tab <- returnTable(t4m3, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                                    "US Pressure", "Rule of Law",
                                    "Rule of Law Nonlinear",
                                    "Ratification"),
                   xtable=TRUE, caption= "Table 4, Model 1.")
tab
sink()

pdf(file = "figures/T4M3.pdf", width = 10, height = 4, family = "Helvetica") 
par(mfrow=c(1,2))
Greg::plotHR(t4m3, xlim=c(0,60), ylim=c(.5, 4), rug="ticks", main="Spatial Spline", xlab="Number of Roads to Criminalized Neighbors")
Greg::plotHR(t4m3,term=3, rug="ticks", xlim=c(-2.5,2.5), xlab="Rule of Law",
             main="Rule of Law")
dev.off()

diagnose(t4m3)


##
#internal
##
form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2)+ uspressure_1 + 
  pspline(rule_of_law, df=3) + rattocprot + cluster(country)
t4m4 <- coxph(form, data=data[data$internal_trafficking==1,], x=TRUE)


sink(file = "tex/T4M4.tex")
tab <- returnTable(t4m4, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                                    "US Pressure", "Rule of Law",
                                    "Rule of Law Nonlinear",
                                    "Ratification"),
                   xtable=TRUE, caption= "Table 4, Model 1.")
tab
sink()

pdf(file = "figures/T4M4.pdf", width = 10, height = 4, family = "Helvetica") 
par(mfrow=c(1,2))
Greg::plotHR(t4m4, xlim=c(0,60), ylim=c(.5, 4), rug="ticks", main="Spatial Spline", xlab="Number of Roads to Criminalized Neighbors")
Greg::plotHR(t4m4,term=3, rug="ticks", xlim=c(-2.5,2.5), xlab="Rule of Law",
             main="Rule of Law")
dev.off()




########
# Figure 7: Weighting neighbors by type.
########

modeld <- coxph(Surv(t0,t,fail) ~ pspline(originSum, df=2) + 
                  pspline(destinationSum, df=2)+ 
                  transitSum +
                  pspline(rule_of_law, df=3) +
                  uspressure_1 + rattocprot  + 
                  cluster(country), data=data)

sink(file = "tex/T5.tex")
tab <- returnTable(modeld, rownames=c("Origin", "Origin Nonlinear",
                                    "Destination", "Destination Nonlinear",
                                    "Transit", 
                                    "Rule of Law", "Rule of Law Nonlinear",
                                    "US Pressure", "Ratification"),
                   xtable=TRUE, caption= "Results Table for Figure 7")
tab
sink()

pdf(file = "figures/Figure7plus.pdf", width = 15, height = 4) 
par(mfrow=c(1,3))
Greg::plotHR(modeld, xlim=c(0,60), ylim=c(.5, 8), rug="ticks", 
             main="Origin Neighbors", 
             xlab="Number of Roads to Criminalized Origin Neighbors")
Greg::plotHR(modeld, term=2, xlim=c(0,60), ylim=c(.5, 8), rug="ticks",
             main="Destination Neighbors",
             xlab="Number of Roads to Criminalized Destination Neighbors")
Greg::plotHR(modeld,term=4, rug="ticks", xlim=c(-2.5,2.5), xlab="Rule of Law",
             main="Rule of Law")
dev.off()

